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TURBOMACHINERY APPLICATIONS 


R.V. Chima and M.-S. Liou 
National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


Abstract 

Many turbomachinery CFD codes use second-order 
central-difference (C-D) schemes with artificial viscos- 
ity to control point decoupling and to capture shocks. 
While C-D schemes generally give accurate results, they 
can also exhibit minor numerical problems including 
overshoots at shocks and at the edges of viscous layers, 
and smearing of shocks and other flow features. In an 
effort to improve predictive capability for turbomachin- 
ery problems, two C-D codes developed by Chima, 
RVCQ3D and Swift, were modified by the addition of 
two upwind schemes: the AUSM + scheme developed by 
Liou, et ah, and the H-CUSP scheme developed by Tat- 
sumi, et al. Details of the C-D scheme and the two 
upwind schemes are described, and results of three test 
cases are shown. Results for a 2-D transonic turbine 
vane showed that the upwind schemes eliminated vis- 
cous layer overshoots. Results for a 3-D turbine vane 
showed that the upwind schemes gave improved predic- 
tions of exit flow angles and losses, although the H- 
CUSP scheme predicted slightly higher losses than the 
other schemes. Results for a 3-D supersonic compressor 
(NASA rotor 37) showed that the AUSM + scheme pre- 
dicted exit distributions of total pressure and tempera- 
ture that are not generally captured by C-D codes. All 
schemes showed similar convergence rates, but the 
upwind schemes required considerably more CPU time 
per iteration. 

Introduction 

Turbomachinery blades are usually designed with 
proprietary design codes and are heavily analyzed with 
computational fluid dynamics (CFD ) codes before com- 
mitting to manufacture. However, turbomachinery 
designers often distrust absolute performance predic- 
tions and rely only on changes in predicted performance 
between designs. This practice suggests that the accu- 
racy of CFD codes can still be improved. 

In 1994 ASME and IGTI sponsored a blind test 
case for turbomachinery CFD codes at the 39th Interna- 
tional Gas Turbine Conference held in The Hague 
(unpublished.) The same test case was later adopted by 


the AGARD Propulsion and Energetics Panel Working 
Group 26 as a test case for examining effects of grid and 
turbulence model on solution accuracy. ’ Sixteen dif- 
ferent CFD codes were used to predict the performance 
of a transonic compressor rotor designated NASA rotor 
37. 3 4 One operating point at 98 percent of maximum 
flow was examined in detail. Predicted pressure ratios 
varied by nearly 10 percent, and predicted efficiencies 
varied by about 6 points. In general, pressure ratios were 
too high and the efficiencies were too low. The large 
variations in results again suggests that the codes can 
still be improved. 

Most of the codes used for these test cases used sec- 
ond-order central-difference (C-D) schemes with artifi- 
cial viscosity to control point decoupling and to capture 
shocks. While C-D schemes generally give accurate 
answers, they can also exhibit some minor numerical 
problems. Shock smearing and overshoots are well 
known and can be minimized by switching off fourth- 
difference dissipation at shocks. Overshoots at the edge 
of viscous layers are less well known but were shown in 
refs. 5 and 6. Many researchers have speculated that 
artificial viscosity may smear out other flow features, 
but this can be difficult to demonstrate. 

Other work has shown that improved artificial vis- 
cosity schemes or upwind schemes can give better accu- 
racy than standard C-D schemes. In ref. 5 Tweedt, 
Chima, and Turkel compared two artificial viscosity 
schemes in a C-D code. The first was a standard artifi- 
cial viscosity scheme with blended second and fourth 
differences and Eigenvalue scaling. The second was 
the Symmetric Limited Positive (SLIP) flux limiter 
which is the low-speed part of the more general Convec- 
tive Upward Split Pressure (CUSP) schemes developed 
by Tatsumi, Martinelli, and Jameson. 10 ' 11 It was shown 
that the SLIP formulation gave better resolution of lami- 
nar boundary layer velocity profiles and better predic- 
tions of performance of a low-speed centrifugal impeller 
than the standard formulation. 

In several papers, Liou and others have developed 
the Advection Upstream Splitting Method (AUSM) 
family of upwind schemes and applied them to many 


NASA/TM— 2003-212457 


1 



aerodynamic problems ranging from 1-D shock tube 
problems to 3-D multi -element wings. These appli- 
cations have shown that the AUSM schemes have excel- 
lent shock-capturing properties and give very accurate 
results for a wide variety of problems. 

In the present work the H-CUSP and AUSM + 
schemes were added to two C-D turbomachinery analy- 
sis codes developed by Chima, RVCQ3D 616 and 
Swift. 17-19 Three turbomachinery blades were analyzed 
and the results were compared to experimental data. In 
each case the upwind schemes gave significant improve- 
ments over the C-D scheme. Results for a 2-D transonic 
turbine vane showed that the upwind schemes elimi- 
nated viscous layer overshoots. Results for a 3-D turbine 
vane showed improvements in predicted exit flow angle 
and loss profiles with the upwind schemes. Finally, 
results for a 3-D supersonic compressor (NASA rotor 
37) showed that the AUSM + scheme gave large 
improvements in the prediction of exit total pressure and 
total temperature profiles. 

CFD Codes 

Swift 

The Swift code is a multiblock Navier-Stokes anal- 
ysis code for turbomachinery blade rows. 17-19 The code 
solves the Navier-Stokes equations on body-fitted grids 
using an explicit finite-difference scheme. It includes 
viscous terms in the blade-to-blade and hub-to-tip direc- 
tions, but neglects them in the streamwise direction 
using the thin-layer approximation. Two turbulence 
models were used: the Baldwin-Lomax model, and 
Wilcox’s k-co model 6 71 with Menter’s shear stress 
transport (SST) modification. 77 

The baseline code used C-D’s for the fluxes, and 
scalar artificial dissipation to capture shocks and to con- 
trol point decoupling. 17 Eigenvalue scaling was used to 
scale the artificial dissipation directionally on stretched 
grids. 8 ' 9 An explicit, four-stage Runge-Kutta scheme 7 
was used to solve the flow equations. To accelerate con- 
vergence to a steady state, all calculations were run at a 
Courant numbers around 5.6 using a spatially- varying 
time step and implicit residual smoothing. The Eigen- 
value scaling was also used to scale the implicit smooth- 
ing coefficients. 

RVCQ3D 

The quasi-three-dimensional turbomachinery anal- 
ysis code RVCQ3D developed by Chima 616 was used to 
develop and test the upwind schemes before attempting 
3-D calculations. RVCQ3D solves the thin-layer Navier- 


Stokes equations on a blade-to-blade plane. Radius 
change, stream surface thickness, and rotation can all be 
modeled. The differencing scheme, artificial dissipation, 
and solution algorithms were all similar to those 
described previously for the Swift code. 

Governing Equations 

The Navier-Stokes equations were written in a Car- 
tesian ( x , y, z) coordinate system rotating about the x- 
axis with angular velocity 12. The equations were trans- 
formed to a curvilinear (£, q, q) system using standard 
techniques, and all viscous terms in the ^-direction were 
dropped using the thin-layer approximation. The result- 
ing equations are: 

r) t q + J[d^E + r) n F + d q G - Re~ l (3,^,, + 3 ? G t ,)] = H (1) 

where q = [p, p«, pv, p w, e] T is the vector of conserva- 
tion variables, and 

P U' 

P uU’ + ^p 

E = J~ l pvU' + t^ y p (2) 

pwt/' + ^p 

p h Q U' + 0.^ eP 

etc. are the inviscid fluxes, F and G v are viscous fluxes, 
and H is a source term due to rotation. In (1) and (2) 

e = p\C v T + i(i i 2 + v 2 + w 2 )J is the total internal energy 

and h 0 = (e + p )/ p is the total enthalpy. The full equa- 
tions are given in ref. 17. 

In equation (2) U' is the contravariant velocity 
component in the relative frame of reference. Using 
primes to denote relative velocities, 

u' = ^ 

v' = v - Slz (3) 

w' = w + £2y 

Rearranging terms gives 

U' = (%u + %V + %w)-Si(%7-%y) 

' „ (4) 

= U-Q.^ 

Metric terms ^ , etc. are evaluated at grid points 

using a conservative, centered scheme. 5 The metric 
terms (including £ e ) are averaged to i ± 1/2 for the 
upwind schemes. The Jacobian term J can usually be 


NASA/TM— 2003-212457 


2 



combined with the metric terms and will be neglected 
here. 

Equation (1) is solved in the following form: 

d,q = -J[R.-(R V + D)] (5) 

where Rj is the inviscid residual, R v is the viscous resid- 
ual, and D is a numerical dissipation operator. 

Artificial Dissipation and Upwind Schemes 
Baseline Scheme 

The baseline numerical scheme uses standard cen- 
tral differences for the flux terms 3 , and a scalar arti- 
ficial dissipation D. D is written as the sum of a second 
and forth difference operator in each direction: 

Dq = (Dp+D +D )q 

9 5 ( 6 ) 
D^q = 3^[C^(£ 2 3^<7 - e 4 d^q)] 


sound, and q s is the inverse of the spacing normal to the 
surface. 

H-CUSP Scheme 

The Symmetric Limited Positive (SLIP) scheme 
was introduced by Tatsumi, Martinelli, and Jameson in 
refs. 10 and 11. The scheme uses flux-limited dissipa- 
tion to produce a non-oscillatory scheme. 

The Convective Upward Split Pressure (CUSP) 
scheme was also introduced in refs. 10 and 11. The 
CUSP scheme was developed as a flux-split scheme 
similar to the AUSM scheme; however, it was imple- 
mented as a dissipative flux added to a C-D flux. For 
computational efficiency the dissipative fluxes can be 
updated less often than the C-D fluxes. The E-CUSP 
formulation bases the dissipative fluxes on the internal 
energy, while the H-CUSP formulation is based on stag- 
nation enthalpy. The H-CUSP formulation was used 
here. 


where £, and e 4 are coefficients given by 


For the H-CUSP scheme, the artificial dissipation is 
written as: 


e 2 = ^2 mfl *(V/_i>v,>v,. + i,v ( . + 2 ) 
e 4 = max( 0, K 4 - e 2 )/ ; 

K 2 = 1/4 

K 4 = (0.25 to 1)/16 
v ( is a pressure sensor for shocks 


= \Pi-i- 2 Pi + Pi + i\ 

'' \p i _ l + 2 Pi + Pi+1 


( 8 ) 


and /. is a ramping function that reduces the dissipation 

linearly with grid index near solid surfaces (typically by 
a factor of 0.05 at the wall) to minimize effects on skin 
friction. 


Experience has shown that more dissipation is often 
needed along the long side of highly-stretched cells. 
Thus equation (6) includes an Eigenvalue scaling coeffi- 
cient Cf , originally proposed by Martinelli, et al. x Here 

a modification proposed by Kunz, et al. 9 was used. 



(9) 


D ifl = (^+i/2- rf ; + i/ 2 )<7 0°) 

where 

d i+ 1/2 = \a*c%,{q R -q L ) + \${E R -E L ) ( 11 ) 

The first term is a difference of the conservation 
variables q. If q R = q j+l and q L = q i , the term 

becomes a first-order artificial dissipation. If q R and q L 

are evaluated using the SLIP limiter described later, the 
term becomes third order in smooth regions of the flow 
and first order near shocks, similar to baseline scheme. 
The second term is a difference of the fluxes E. It is 
added to give a true upwind scheme for supersonic flow. 
Switching terms a* and p are devised to use the first 
term for low speeds and the second term for M > 1 , with 
a continuous blending in between. 

Liou and Steffen proposed a decomposition of the 
flux E into a convective term and a pressure term. An 
equivalent splitting is used for the CUSP schemes. 

E = U'y + pg 
V = [p, p«, pv, pw, ph Q \ T 

g = 


In (9,) /= is the maximum Eigenvalue (i.e., the spectral Using this decomposition, (11) becomes: 
radius) of the inviscid flux Jacobian, c is the speed of 


NASA/TM— 2003-212457 


3 



d i + i/i = \**c%te R -<h) ( 13 ^ 

+ 2 p[(^R U L ¥l) + Sj+ 1/2 (Pr~Pi)1 
An interface relative Mach number is defined by: 


M = ^L- (14) 

ct, s 

where the tilde indicates Roe averaging. 


u 


c 


JP~L U L + JPr 11 R 

— — — , etc., and 

JPl + JPr 

J(Y- 1 )(h 0 -^(u 2 + v 2 + w 2 )j 


Then the switching function p is given by: 


(15) 


P = 


max{ 0, 2M - 1 ) for 0 < M < 1 
min(0,2M + 1) for -1<M<0 
sign(M) for \M\ > 1 


(16) 


which can be coded conveniently as 
fit = 2|Af| - 1 

(17) 

P = signal, M) x min[l, max(Q, f M )] 

Tatsumi, et al. showed that one -point shocks could 
be obtained if switching function a*c is given by: 

a*c = ac- pt/' (18) 

where 


a = r + max(\M\, r'M 0 ) 
r + = max[l, (T— 2|M|)r^] 
r = min(l,r 
r t, = + 


(19) 


M q is a cutoff Mach number, typically taken as 
M q ~ 0.1 x min( 1, M max ) , where M is the largest rela- 
tive Mach number expected in the flow field. Swanson, 
et al. showed that the CUSP schemes also benefit from 
increased dissipation along the long side of stretched 
cells. 23 They suggested Eigenvalue scaling terms 
r + and r . Here r + has been modified by the addition of 
the (1 -2|M|) term. Coefficients a and P are shown in 
figure 1. 



Figure 1 — Functions a and p for the H-CUSP 
scheme 

AUSM + Scheme 

The Advection Upstream Splitting Method 
(AUSM) scheme was introduced by Liou and Steffen in 
1991 12 -p^g AIJSM scheme defines a cell interface 
Mach number based on characteristic speeds from the 
neighboring cells. The interface Mach number is used to 
determine the upwind extrapolation for the convective 
part of the inviscid fluxes. A separate splitting is used 
for the pressure terms. Generalized Mach number and 
pressure splitting functions were described by Liou 13 
and the new scheme was termed ASUM + . The AUSM + 
scheme was shown to have several desirable properties: 
1, it gives exact resolution of 1-D contact and shock dis- 
continuities, 2, it preserves positivity of scalar quanti- 
ties, and 3, it is free of oscillations at stationary and 
moving shocks. 

The AUSM + scheme avoids an explicit artificial 
dissipation, and differences the fluxes directly using: 

= E !+ 1/2 -Ej_ i/2 (20) 

A flux decomposition similar to (12) is used to 
write 


E = p C/'4> + pg 

( 21 ) 

0 = \|//p = [1 ,u,v,w,h 0 Y 

Here p U' is the mass flux across a cell interface. It can 
be written as: 

m = p U' = p c% s ¥r = P c\ s M (22) 

c ^s 

where c = {c R + c L )/ 2 is the average speed of sound, 
and M is the relative interface Mach number. 


The fluxes are differenced using: 
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>0 


_ J PA M i + 1/2 + 8 Pi + 1/2 if M i + 1/2 
[p s c^ s M ( + 1/2*1)^ + g// (+ 1/2 else (23) 


where and g are evaluated at i + 1 /2 . Note that the 
average speed of sound c has been replaced with a 
numerical speed of sound c which is described later. 

The interface Mach number and pressure are evalu- 
ated using weighted averages of the left and right states. 
Defining left and right Mach numbers based on c as: 




UL' 

^ s' L.R 


then M j + 1/2 and p i+l/2 are given by: 


(24) 


M /+ i/ 2 = M + M + D p 
Pi + 1/2 = P + Pl + p 'Pr + D v 


(25) 


M ± and P ± are functions of M, and M 


if \M l \ 

< 1 

M + 

= m 2l , 

P + =P 5 L 

else 


M + 

= m il , 

sT 

ii 

+ 

if \M r \ 

< 1 

M " 

= M 2 R ’ 

P =P 5 R 

else 


M 

= M lR , 

P' = M 1R . 


(26) 


M 


2(7 R) are second order polynomials in M 


l.r ■ 


M 1L = \(M l+ 1)2 


M 


2 R 


\(M r - l) 2 , (27) 


M 


1(i R) are directional switching functions: 


M IL = max(M L ,0) R = min(M R ,Q), 


(28) 


and P 


5 (L.R) 


are fifth order polynomials in M 


L.R ■ 


P 5L = M 2L (2-M L ) + ±M L (Ml-l)* 

(29) 

P 5R = -M 2S (2 + M S )^^M S (M|- 1) 2 

Plots of M~ and P ± are shown in figure 2. 

The AUSM schemes define the mass flux across a 
cell interface in terms of split Mach numbers and a com- 
mon interface speed of sound, c in equation (22.) How- 
ever, Liou and Edwards showed that the interface speed 
of sound c in equation (24) could be chosen arbitrarily 
without affecting the shock-capturing properties of the 



Figure 2 — Functions M ± and /^for the AUSM + 
scheme 


scheme. 14 They proposed using a “numerical speed of 
sound” that effectively scales the numerical dissipation 
with the local flow speed |i/| instead of the local sound 
speed c as M — > 0 . In other words, the numerical speed 
of sound goes to zero with the local Mach number. They 
showed that the numerical speed of sound gave appro- 
priate amounts of dissipation, even when used with pre- 
conditioning methods at very low speeds. 

The numerical speed of sound is given by: 


c = fc 


f = 


+ 4M} 

1 +Ml 


M = l -(M L + M R ) 

M t . = min[l, max(\M\, M q )] 


(30) 


where f is a scaling factor, M is an average interface 
Mach number, and M, is the local relative Mach num- 
ber limited between a cutoff Mach number M Q and 1 . 
M 0 is typically taken as (0.2 to 0.5) x min(M max , 1) . 

In equations (25,) D p and D v are diffusive terms 
that have been introduced to ensure pressure-velocity 
coupling at low speeds. D p is a pressure-diffusion term 

that was introduced by Liou and Edwards. 14 The term 
was originally added to the mass flux, but here it was 
recast as a modification to M j+l/ 2 , and all density 

terms were cancelled. The term was also reduced by a 
factor of two by numerical experimentation. The result 
is: 
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D = 1 A m (P L ~Pr) 

P 4 mU Pl + Pr ) (31) 

AM = (M 2L -M lL )-(M 2R -M lR ) 

Finally D v is a velocity-diffusion term that was 
introduced by Liou. 15 It is given by: 

D v = -P + p'(^1^~cHm r - Ml ) (32) 

Limiters 


SST k-co Turbulence Model 

Results for a transonic compressor rotor shown 
later used the SST k-co turbulence model. Wilcox’s 
baseline k-co model was described in ref. 21, and the 
implementation of the model in RVCQ3D was described 
in ref. 6. The shear stress transport (SST) model was 
developed by Menter in ref. 22, and is described below. 

The SST model is based on Bradshaw’s assumption 
that the shear stress in a boundary layer is proportional 
to k. Menter showed that this could be added to the base- 
line model as a modification to the turbulent viscosity. 


SLIP Limiter 

For the H-CUSP scheme the right and left states 
were calculated using the SLIP limiter. 10 For left and 
right states a and b, the limiter is defined by: 

R(a,b) = 1- a ~ b 2 , (33) 

M + \b\ 

and a limited average is defined by 

L(a,b) = R(a,b){^y (34) 

The conservation variables q are interpolated to 
i ± 1/2 using: 

c 1l = ( li + A< ?! + 3/ 2 ) 

c 1r = < li+ 1 + 2^ A< 7;-l/2> A< 7; + 3/2) 
a< 7;- t/2 = C U~ C U - 1 

van Albada Limiter 

For the AUSM + scheme the right and left states 
were calculated using the van Albada limiter.” The lim- 
iter is defined by: 


S( a b ) = ( a2 + e ) v + ( fc2 + £ ) t ' 
( a 2 + e) + ( b 2 + e) 


(36) 


where e = 10 7 . The primitive variables 
co = [j^ Vj w , ; p ] T are interpolated to i± 1/2 using: 


"/. = co i +is , (Aco,._ 1 /2 ) Aco,. + i/ 2 ) 

a R = ®i+ 1 ~ 2 S ^ (0 '- 1/2 ’ ACO i'+ 1/2 ^ 


(37) 


a^k 

V? max(a^(0, SIF 0 ) ^ 

where a t = 0.3 , and O is the magnitude of the vortic- 

ity. The first term in the denominator recovers the base- 
line model and the second term gives the SST model. 
Since Bradshaw's assumption does not necessarily hold 
in free-shear layers, F 2 is a blending function that turns 
the SST model off away from the wall. 


F 2 = tanh(arg?) 


arg 2 = max 


2 Jk 400 v 3 

0.09 coy coy 2 J 


(39) 


Menter has shown that the SST model gives excel- 
lent results for adverse pressure gradients. The one dis- 
advantage to the model is that it requires the distance to 
the wall y. 

Results 


2-D Transonic Turbine Vane 

A transonic turbine vane tested by Arts, et al. was 
computed as a 2-D test case. The vane was tested exper- 
imentally in the Isentropic Light Piston Compression 
Tube Facility at the von Karman Institute. The facility 
has independent control over the exit Reynolds number, 
the exit isentropic Mach number, M, ■ , and the inlet 

turbulence intensity. Surface pressures were measured 
with static taps, and wake total pressure profiles were 
measured with a high-speed traversing probe. 

For the computations a C-type grid was used with 

383 x49 (C, y) points. The grid spacing gave y + < 1.5 
over most of the vane. The grid size was found to give 
good resolution of the suction surface shock and surface 
heat transfer in ref. 6. Solutions were run for a case with 
an exit Mach number of M 2 is = 1.0 using the C-D, H- 

CUSP, and the AUSM + schemes and the Baldwin- 
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Lomax turbulence model. All solutions were run using a 
4-stage Runge-Kutta scheme at a Courant number of 
5.6. Dissipative terms were evaluated after the first 
stage, and the turbulence model was updated every five 
iterations. Convergence rates were similar for all 
schemes, with mass flow error and total pressure loss 
converged to 0.3 percent or better in 2000 iterations. On 
an SGI Octane workstation the C-D solution took 188 
sec. The H-CUSP solution took 1.27 times longer, and 
the AUSM + solution took 2.65 times longer than the C- 
D scheme. 

Computed Mach contours are shown in fig. 3. The 
heavy black line is M = 1.0 and the contour increment 
is 0.05. The flow accelerates from M = 0.15 at the inlet 
to M = 1.2 on the suction surface. A normal shock 
reduces the Mach number to M~ 1.0 at the exit. 
Enlargements of the shock, trailing edge, and wake 
computed with the three different schemes are shown in 
fig. 4. The C-D results show some oscillations around 
the shock and more severe oscillations around the wake. 
The H-CUSP scheme eliminates most of the oscilla- 
tions, although some are visible in the core flow. The 
AUSM + results show a very clean shock and are com- 
pletely non-oscillatory. 

Computed distributions of isentropic surface Mach 
number are compared to experimental data in fig. 5. ' All 
schemes agree very well with the experimental data. 

Computed wake profiles located 43 percent of axial 
chord downstream of the trailing edge are compared to 
the experimental data (digitized manually from ref. 25) 
in fig. 6. The C-D results show the same oscillations in 
total pressure at the edge of the wake that were seen in 
the contour plots in fig. 4. Neither of the upwind 
schemes shows oscillations. The AUSM + results agree 
very well with the experimental data, but the H-CUSP 
results show slightly too much wake decay and free- 
stream loss. 

3-D Subsonic Turbine Vane 

An annular turbine vane that was tested experimen- 
tally by Goldman and McLallin at NASA Glenn 
Research Center 26 was used as a 3-D turbine test case. 
A C-type computational grid was used, with 
97x32x33 (C, 0, r) points. The grid spacing gave 

y + = 0(5) over most of the vane. Although the grid 
was rather coarse, it gave reasonably accurate predic- 

f. In all figures the C-D results are labeled “Baseline” 
and are shown with a solid black line, AUSM + results 
are shown with a dashed red line, and H-CUSP results 
are shown with a dotted blue line. 


tions of vane performance with quick turnaround. Solu- 
tions were run using the C-D, H-CUSP, and the AUSM + 
schemes and the Baldwin-Lomax turbulence model. The 
iterative scheme was the same as that used for the 2-D 
case. Convergence rates were similar for all schemes, 
with the maximum residual reduced about four orders of 
magnitude in 1500 iterations. Total pressure losses were 
converged to four digits. 

The solutions were run on the Cray SVlex com- 
puter at NASA Ames Research Center (Bright.) The C- 
D solution took about 1/2 hour, or about six minutes of 
wall clock time using six processors. The H-CUSP 
scheme took 1.47 times longer than the C-D scheme, 
and the AUSM + scheme took 1.57 times longer than the 
C-D scheme. 

Computed pressure contours on the blade surfaces 
are shown in fig. 7. The blade profile is uniform along 
the span, and the pressure distribution is nearly uniform. 
The flow accelerates from M = 0.21 at the inlet to 
M = 0.73 at the exit. 

Figure 8 compares measured and calculated con- 
tours of kinetic energy efficiency across the wake at a 
distance of 1/3 axial chord downstream of the trailing 
edge. The kinetic energy efficiency is defined by: 


where Q is the velocity, Tq is the total temperature, and 
T is the static temperature. The C-D scheme smears 
many of the details of the wake. The H-CUSP scheme 
captures the wake shape better, showing underturned, 
high loss regions near the endwalls due to secondary 
flows. The AUSM + scheme overexaggerates the wake 
shape; however, subsequent results will show that the 
AUSM + scheme gives the best quantitative agreement 
with experiment. 

The spanwise variation of mixed out total pressure 
loss coefficient 1 - Po^/Pq,-,, downstream of the vanes 

is shown in fig. 9. The C-D results show little detail 
along the span. The H-CUSP results show some detail 
near the tip but too much loss near the hub. The AUSM + 
results show good qualitative agreement with the data 
along the entire span. All results show higher losses than 
the data at midspan. The midspan loss does not improve 
with increasing grid resolution, and may be due to poor 
modeling of the round trailing edge. 

The spanwise variation of flow angle downstream 
of the vanes is shown in fig. 10. The C-D results show 
nearly uniform flow angle along the span, and the H- 
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CUSP results are only slightly better. The AUSM + 
results show excellent agreement with the data along the 
entire span. 

3-D Compressor Rotor 

A low aspect ratio transonic inlet rotor for a core 
compressor, designated NASA rotor 37, was used as a 3- 
D compressor test case. The rotor was originally 
designed and tested at NASA Glenn Research Center in 
the late 1970’s by Reid and Moore. 3 ' 4 It has 36 multiple- 
circular-arc blades and a design pressure ratio of 2.106 
at a mass flow of 20.19 kg/sec. 

The rotor was re-tested in a single-stage compressor 
facility at NASA Glenn. The test facility was described 

97 9R 

by Suder, et al. ’ Radial distributions of static and 
total pressure, total temperature, and flow angle were 
measured at two axial stations located 4.19 cm upstream 
and 10.19 cm downstream of the blade hub leading 
edge. 

These measurements were used for the ASME/IGTI 
blind test case and the AGARD test case for turboma- 
chinery CFD codes. Calculations from sixteen differ- 
ent CFD codes were compared to the measurements. 
Two details of the measurements proved to be difficult 
to predict: First, most codes overpredicted the overall 
pressure and temperature ratios, and underpredicted the 
efficiency. Second, most codes failed to predict the 
radial distributions of P Q and T 0 downstream of the 

rotor. Measured distributions show deficits in these 
quantities near 20 percent span, but most codes showed 
fairly linear radial distributions. 

Two researchers predicted these distributions cor- 
rectly: Hah using his HAH3D code with a pressure- 
based, high-order upwind difference scheme and a k-e 
turbulence model," and Weber using the OVERFLOW 
code with a Roe upwind scheme and the Spalart-Alma- 
ras turbulence model. 1,2 Most of the other codes used 
for the test case used C-D schemes with artificial viscos- 
ity and a variety of turbulence models. 

Hah believed that the deficits in total conditions 
were due to a corner stall. 29 Alternatively, Shabbir, et al. 
proposed that flow leakage between the centerbody and 
rotor disk could generate enough blockage to produce 
the deficits. 30 In this paper we suggest that the C-D 
schemes used in most codes smear details of the P 0 and 

T q distributions, while the upwind schemes used previ- 
ously by Hah and Weber, and now in this work provide 
increased accuracy that gives better agreement with the 
experimental data. 


A multiblock grid was used for the present calcula- 
tions (fig. 11.) An H-type grid was used upstream of the 
blade with 45 x 34 x 63 (x, 6, /•) points. A periodic C- 
type grid was used around the blade with 259 x 46 x 63 
points. The grid spacing at the blade and endwalls was 

4xlCT 4 cm, giving y + = 2 - 4 at the surfaces. The blade- 
to-blade grid was optimized in a grid refinement study per- 
formed for the ASME/IGTI blind test case (unpublished.) 
The inlet and exit of the grid were coincident with the 
measurement stations described earlier. An O-type grid 
was used above the tip of the blade with 199 x 13 x 13 
points (13 points across the gap.) The total grid had 
869,011 points, which is 3-4 times finer than the grids 
recommended by Dunham, et al. 

C-D/Baldwin-Lomax calculations were run previ- 
ously for the ASME blind test case. 18 Some of these 
results are included later for comparison. 

The present results were computed using the SST k- 
ft) turbulence model. Preliminary results using the base- 
line k-co model showed that the AUSM/k-co scheme pre- 
dicted higher pressure ratios than the AUSM/Baldwin- 
Lomax scheme. The reason was unclear, but seemed to 
be related to better resolution of the shock/boundary 
layer interaction on the casing. Menter’s SST k-co model 
was then added to the baseline k-co model. Pressure 
ratios predicted with the AUSM/SST k-co scheme 
agreed closely with the AUSM/Baldwin-Lomax results 
and are presented here. 

The ASUM + scheme was used to calculate several 
operating points. The C-D and H-CUSP schemes were 
each used to compute one operating point at 98.7 per- 
cent max flow. All calculations were run with a four- 
stage Runge-Kutta scheme at a Courant number of 5.5. 
Artificial and physical dissipation terms were evaluated 
at stages 1 and 2. The turbulence model was updated 
every two iterations. The calculations were typically run 
3, 000 iterations to ensure convergence of the mass flow 
error and total pressure ratio to about 0.01 percent. The 
total CPU time on the Cray SVlex computer was about 
10 hours per case for the C-D scheme, but on six proces- 
sors the wall clock time was roughly 1.8 hours. The H- 
CUSP scheme took 1.20 times longer than the C-D 
scheme, and the AUSM + scheme took 1.24 times longer 
than the C-D scheme. 

Figure 12 shows computed contours of relative 
Mach number at 73 percent span at 98.7 percent max 
flow. The heavy black contour is M = 1 and the contour 
increment is 0.05. An oblique shock system runs 
upstream of the blade and across the passage, where it 
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merges with a normal shock. AUSM + results are shown, 
but the C-D and H-CUSP results look similar. 

Computed maps of total pressure ratio, total tem- 
perature ratio, and adiabatic efficiency versus mass flow 
are shown in fig. 13. The dotted black line shows the C- 
D/Baldwin-Lomax results reported in ref. 18, with the 
one CD/SST k-co result added. The blue triangles show 
the H-CUSP solution, and the red triangles show the 
AUSM + solutions. No attempt was made to determine 
the numerical stall point with the AUSM + scheme. All 
schemes overpredict the pressure and temperature 
ratios, but give very good predictions of the adiabatic 
efficiency. 

Figures 14 - 16 compare radial profiles of total 
pressure, total temperature, and adiabatic efficiency 
downstream of the rotor with experimental data taken at 
98 percent of the maximum flow rate. Hah, et al. showed 
that these profiles were very sensitive to the flow rate, 29 
and that much better agreement was obtained by com- 
paring calculations at about 99 percent flow. The solu- 
tions shown here are all at a flow rate around 98.7 
percent max flow. 

Total pressure profiles are shown in fig. 14. The 
data shows the deficit in P Q below 30 percent span that 

most codes in the ASME/AGARD test case were unable 
to predict. Here the H-CUSP results show a nearly linear 
distribution of P Q along the span that still fits the data 
well overall. The baseline C-D results are similar near 
the tip but show an overshoot near the hub. The AUSM + 
results match the data very well except for a slight over- 
shoot at the hub. Many of the codes in the ASME/ 
AGARD test case showed similar overshoots near the 
hub. 

Total temperature profiles are shown in fig. 15. The 
C-D results are smooth along the span and do not match 
the shape of the measured profile very well. The H- 
CUSP results are similar, but give slightly better resolu- 
tion of the profile shape. The AUSM + results agree very 
closely with the data between 15 and 85 percent span. 
The three schemes give minor differences in predicted 
T q near the hub that are consistent with the overshoots 

in P Q noted above. All three schemes overpredict T Q 

near the tip, which accounts for the high overall temper- 
ature (and pressure) ratios in fig. 13. Almost every code 
in the ASME/AGARD test case also overpredicted T 0 
near the tip, and the reason remains unknown. 

Adiabatic efficiency profiles are shown in fig. 16. 
Here all three schemes give remarkably similar results 


that agree very well with the data below 85 percent 
span. This indicates that loss levels are being predicted 
correctly by the SST k-co turbulence model, except per- 
haps near the casing. 

Conclusions 

Two centrally-differenced (C-D) turbomachinery 
analysis codes developed by Chima, RVCQ3D and 
Swift, were modified by the addition of two upwind 
schemes: the AUSM + scheme developed by Liou, et al. 
and the H-CUSP scheme developed by Tatsumi, et al. 
Several test cases were run to evaluate the effects of the 
differencing schemes on turbomachinery flow predic- 
tions. The upwind schemes gave improvements in the 
predictions over the C-D scheme for every case investi- 
gated. The following results were noted: 

1 . The C-D scheme produced overshoots at the edge of 
viscous layers. These were eliminated by both the 
AUSM + and H-CUSP schemes. 

2. Although the AUSM + and H-CUSP schemes have 
excellent shock capturing properties for model prob- 
lems, all schemes gave comparable shock resolution 
on general grids. 

3. The H-CUSP scheme usually predicted slightly 
lower total pressures (higher losses) than the other 
schemes. 

4. There was no significant difference in convergence 
rates for the three schemes. 

5. The C-D scheme has the lowest operation count and 
required the least CPU time of the three schemes. 
The H-CUSP scheme uses the same inviscid fluxes 
as the C-D scheme but has more complicated dissi- 
pative fluxes and is therefore slower. In both 
schemes the dissipative fluxes can be updated after 
the first one or two stages of a multistage Runge- 
Kutta scheme to save time. The AUSM + scheme has 
the highest operation count and was updated every 
stage, so it was the slowest of the three schemes. For 
a 2-D problem the H-CUSP scheme was 1.27 times 
slower than the C-D scheme and the AUSM + 
scheme was 2.6 times slower. For 3-D problems the 
viscous fluxes and turbulence models require dis- 
proportionately more time than the inviscid fluxes, 
so the AUSM + scheme requires relatively less of the 
overall time. For 3-D problems the H-CUSP scheme 
was 1.20 - 1.47 times slower than the C-D scheme 
and the AUSM + scheme was 1.24 - 1.57 times 
slower. 

6. For a subsonic turbine vane the AUSM + and H- 
CUSP schemes predicted the 3-D wake shape better 
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than the C-D scheme. The AUSM + scheme gave the 
best overall predictions of turning and loss distribu- 
tions. 

7. For a transonic compressor rotor the AUSM + scheme 
predicted deficits in total pressure and total tempera- 
ture that were measured experimentally but were not 
generally predicted by the C-D codes used for the 
ASME/AGARD test case. This result was consistent 
with predictions by Hah and Weber using two other 
upwind codes. We believe that the measured deficits 
in total pressure and total temperature are an intrin- 
sic feature of this rotor blade and not a result of hub 
leakage as suggested by Shabbir, et al. Furthermore, 
we believe that C-D schemes tend to smear out these 
details due to relatively coarse spanwise grids, but 
that upwind schemes are able to capture them prop- 
erly. 
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Figure 3 — Computed Mach contours for the VKI 
turbine vane, AUSM + scheme 



Central-difference scheme 



AUSM + 

Figure 4 — Computed Mach contours at the trailing 
edge using three differencing schemes 



Figure 5 — Computed and measured distributions 
of isentropic Mach number on the vane surface 
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Figure 6 — Computed and measured total pres- 
sure profiles 0.43 chords behind the vane 
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Figure 7 — Computed pressure contours on the 
Goldman turbine vane, AUSM + scheme 



Figure 9 — Computed and measured profiles of P 0 
loss coefficient downstream of the vane 



Figure 10 — Computed and measured profiles of 
flow angle downstream of the vane 



Experimental measurements 
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Figure 8 — Measured and computed profiles of 
kinetic energy efficiency in the vane wake 
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Figure 1 1 — Computational grid for NASA rotor 37 



Figure 1 2 — Computed Mach contours at 73 per- 
cent span, AUSM + scheme 
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Figure 13 — Measured and computed operating 
maps of P 0 , T 0 , and r| for rotor 37 
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Figure 14 — Measured and computed profiles of 
total pressure downstream of the rotor 
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Figure 15 — Measured and computed profiles of 
total temperature downstream of the rotor 
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Figure 16 — Measured and computed profiles of 
adiabatic efficiency downstream of the rotor 
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